Periodic Droplet Formation in Chemically Patterned Microchannels 
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Simulations show that when a phase-separated binary AB fluid is driven to flow past chemically 
patterned substrates in a microchannel, the fluid exhibits unique morphological instabilities. For the 
pattern studied, these instabilities give rise to the simultaneous, periodic formation of monodisperse 
droplets of A-in-B and B-in-A. The system bifurcates between time-independent behavior and dif- 
ferent types of regular, non-decaying oscillations in the structural characteristics. The surprisingly 
complex behavior is observed even in the absence of hydrodynamic interactions and arises from the 
interplay between the fluid flow and patterned substrate, which introduces non-linearity into the 
dynamical system. 



Hydrophilic-hydrophobic patterning is used by a vari- 
ety of biosystems to direct the motion of fluids at sur- 
faces. For example, hydrophobically-hydrophilically pat- 
terned backs help desert beetles to capture water, and 
hydrophobic patches control water permeation in leaves 
U . This motif is also used to steer the motion and reac- 
tion of fluid droplets in liquid microchips j2J and is being 
utilized to design self-cleaning substrates |lj . Despite the 
utility of these designs, there have been surprisingly few 
theoretical studies into the dynamics of fluid flow over 
chemically patterned surfaces. In this study, we examine 
a conceptually simple system where two partially mis- 
cible fluids, A and B, are mechanically driven (by an 
imposed shear) to flow past patterned surfaces within a 
microchannel (see Fig. 1). The system exhibits two dis- 
tinct steady-states; however, in the transition between 
the two states, we uncover intricately complicated be- 
havior, where monodisperse droplets of both A-in-B and 
B-in-A are formed periodically in time (as shown in Fig. 
2), and the confined liquid displays regular, non-decaying 
oscillations in its structural characteristics. Furthermore, 
we isolate points where this system bifurcates between 
time-independent behavior and different types of oscil- 
latory patterns. What is striking is that the observed 
phenomena occur even in the absence of hydrodynamic 
interaction; this is distinct from well-known instabilities 
in fluids 0,0 ■ Given that the system is relatively sim- 
ple, the results suggest that complex transitions between 
well-defined steady-states may well be evident in a broad 
variety of dynamical systems. 

The observed complex oscillatory patterns arise from 
a competition between advection and thermodynamics 
as an imposed Poiseuille flow drives the phase-separated 
fluids to flow over the chemically patterned substrates. 
As shown in Fig. 1, the top and bottom of the mi- 
crochannel are decorated with a checkerboard pattern. 
Each checkerboard is composed of two A(B)-like patches, 
which are preferentially wetted by the A(B) fluid. The 
first B patch (in yellow) is placed in the way of the A 
stream (in blue) and correspondingly, the first A patch 
is located in the path of the B fluid. 
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FIG. 1: Schematic of system 



The binary fluid is characterized by the order param- 
eter <p(r,t) — pyi(r,t) ~ ps(r,t), where Pi(r,t) repre- 
sents the local number density of the i-th component, 
i = A, B. The thermodynamic behavior of the system 
is governed by the coarse-grained free energy functional, 
F = Fq + ^fg, where Fq is the Ginzburg-Landau free 
energy for a binary mixture 
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and a and b are positive constants. We consider the 
fluid to be in the two-phase coexistence regime where 
the equilibrium order parameter for the A(B) phase is 



represents the cost 



ifA/B = iy/a/b. The term £ 
of order parameter gradients. The free energy \T/ S de- 
scribes the interaction of a fluid element at a point f 
with the patterned substrate. Specifically, we take 



<F S = I df I ds[ \v{£) ■ e 



\r — s\ lr a 



(<p(f) - ip(s))' 



(2) 

where the inner integral represents integration over the 
substrates. V(s) = V = const on the patterns and is 
zero otherwise [(j , and tq represents the range of the sub- 
strate potential. We choose <p(s) = ^Pa(b) to introduce 
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A{B) -wetted patches at specific regions of the substrate. 
Through eq (0), the free energy F is reduced when the 
fluid is A(B)-rich near A(B)-like patches. The evolu- 



(a) 
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FIG. 2: Periodic droplet formation. Blue represents A-rich 
and yellow represents B-rich domains. Panels a-c show the 
"front" and panels d-f show the "back" of the channel, with 
droplets of B-in-A and A-in-B, respectively, a, d) t — 31200, 
b, e) t = 32000, c,f) t = 33600. The other parameters are: 
H = 3 ■ 10~ 4 , I = 60, h = 40, r = 5, V = 0.003 

tion of the order parameter for this system is described 
by the Cahn-Hillard equation, in which the flux of tp is 
proportional to the gradient of the chemical potential, 
J v = M V/i, where fj, = ^ an d M is the mobility of the 
order parameter. The imposed Poiseuille flow advects 
variations in ip along the microchannel. In dimensionless 
units 
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V 2 



(3) 



where the length scale is chosen to be equal to the thick- 
ness of the interface between A and B fluids, £i nt = 
y/k/a, and the scale of time is the diffusion time through 
that interface, r = ^i nt 2 jaM 0. The velocity field v 
obeys the Navier-Stokes equation in the over-damped 
limit, which is appropriate for low Reynolds number flow, 
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where p is a Lagrange multiplier that guarantees the in- 
compressibility condition, V • v — 0, and H x = (Pj n — 
Pout)iintTrj / L is the dimensionless form of the imposed 
pressure drop (Pj„ — P ou t) along the channel of length 
L. Because the pressure gradient is applied along the x- 
axis, only the x-component of the vector H is non-zero, 
H x = H . The last term in eq. 0] is the non-dissipative 
part of the stress tensor || (for the above dimensionless 
form see 0); this term represents hydrodynamic inter- 
actions. The constant C — a ■ S,i n t/{a ■ r\ ■ M) depends on 
the fluid properties, such as the shear viscosity, r), inter- 
facial tension, a w kip 2 q /^i nt , and diffusivity aM. The 



value of C determines the importance of hydrodynamic 
interactions for the specific fluid. For a fluid with a high 
viscosity, where C « 1, hydrodynamic interactions can 
be neglected. In this work, we set C = 0; therefore, 
the velocity profile in our system is determined by the 
imposed pressure gradient. Thus, advection in a shear 
flow and diffusion of the fluids to the more wettable A or 
B domains control the evolution of the order parameter 
in the system and are responsible for the observed rich 
behavior ||. 

Equations Q~© ar e discretized and solved nu- 
merically by a cell dynamic system method on 
a 120x40x40 grid. The following boundary con- 
ditions on the walls of the channel are imposed: 
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and 

wall 

k~ x J dsi [V(sj) (ip(r) — <p{^i))]\-_ t - The last condition 
arises explicitly from the minimization of free energy in 
the presence of the substrate potential. At the entry of 
the channel, we have two-stream flow; at the exit, we 

assume free draining flow, i.e., =0. For the ve- 

° ' ' ox T 

locity field, we assume no-slip boundary condition on the 
walls 0. 

We consider the order parameter evolution in the chan- 
nel at different values of H . At low H (low velocities), 
local thermodynamics dominate, and the fluid just mim- 
ics the underlying checkerboard pattern, with small dis- 
tortions in ip caused by the imposed flow. Higher veloc- 
ities lead to more dramatic changes in ip and yield the 
complex interfaces between the A/B fluids shown in Fig- 
ures 2 and 3. This behavior occurs when the scale of 
spatial distortions in ip within the center of the channel 
(where the Poiseuille flow exhibits the maximum veloc- 
ity, t> max ) is comparable to the length of the patch, i.e., 
when v max t l di j f ~ I, where ^i// ~ ^ ^ s ^ ne characteris- 
tic diffusion time over the patch length I. This estimate 
yields a value of H « 10~ 4 . For higher values of H, the 
fluid flows through the middle of the channel with essen- 
tially no distortion, while near the top and the bottom 
substrates, ip(r, t) is governed by the patches' wetting 
properties (see Fig. 3e). 

For the intermediate H values, a competition between 
the preferential wetting interactions and the imposed flow 
leads to fascinating behavior. On one hand, the wetting 
effects cause the fluids to diffuse to the respective patches 
to minimize the free energy and there is the general ten- 
dency to minimize interfacial regions between the A and 
B phases. On the other hand, the imposed flow carries 
fluid away from the favorable patches. If both substrates 
contained just the first half of the checkerboard (a single 
A and B patch), these patches would simply "switch" the 
location of the fluids, and the imposed flow would move 
the switched fluids along the channel. The presence of 
the second set of patches interrupts this flow because 
both the A and B streams again confront incompatible 
domains. In three dimensions, each component can avoid 
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the second unfavorable region by diffusing into the bulk. 
However, recall that the first yellow (B) patch is on the 
blue (A) side (see Fig. 2). Thus, the B fluid can extend 
only so far into the incompatible domain (similar argu- 
ments hold for the A fluid). Each fluid forms "arms" that 
reach from the top and bottom of the walls; these arms 
can join and pinch off to a form a bubble. Our coarse- 
grained modeling allows such a topology change without 
any ad hoc rules. Figures 2a-c show the order parame- 
ter distribution at the "front" of the microchannel; the 
same behavior occurs for the A fluid at the "back" of the 
channel (Figs. 2d-f). 




FIG. 3: Order parameter distribution at steady state for 
different H: a) H = 2.6 • fO" 4 , b) H = 2.7 • 10" 4 , c) 
H = 2.79 ■ 1(T 4 , d) H = 3.2 ■ 1(T 4 , e) H = 3.44 ■ 1(T 4 . 
All other parameters are the same as in Fig 2. 

Between the limiting cases in Figures 3a and e for 
relatively low and high H, respectively, three different 
types of behavior have been observed: two types of pe- 
riodic behavior and a time-independent state asymmet- 
ric with respect to the top and bottom substrates (Fig. 
3b); which of the two morphologies shown in Fig. 3b 
is actually realized for fixed H depends on the noise in 
the system JfJ. The periodic cases exhibit "symmetric" 
(Fig. 3c) and "asymmetric" oscillations (Fig. 3d). The 
arms in the figures and all the periodic behavior develop 
mainly near the sidewalls. In the middle of the box, the 
Pouseuille velocity field has a maximum and advection 
prevails. Near the wall, however, the velocity is much 
smaller and diffusion dominates, allowing the arms to 
move upward (downward) and join. 

To analyze the complex dynamics, we define a param- 
eter that characterizes the integrated changes in ip(r,t) 
near the sidewalls: 



Bi(t) = yj t) ~ <p(f, 0))| , (5) 

here i = top, hot indicates whether we integrate over the 
top (see red dashed box in Fig. I) or bottom half of the 
sidewall region of volume T^.(We choose the thickness of 
this region as h/8.) 




t, time 

FIG. 4: Evolution of Bt op and Bbot for cases shown in Fig 3, 
a, c, d and e 



The evolution of Bt op (t) and Bbot(t) for different val- 
ues of H is shown in Fig. 4. Note that Bt op (t) = Bbot(t) 
for all the symmetric cases (Figs. 3 a, c and e). A case of 
asymmetric oscillations, where B top (t) ^ Bb ot {t), is plot- 
ted in red. The maxima in the curves correspond to the 
largest distortions (where the bubbles are biggest); the 
minima correspond to the structure where the arms are 
separated by the greatest distance. The two curves for 
the asymmetric case (Fig. 3d) are similar to each other, 
but there is a phase shift between them. Each of these 
red curves displays two maxima and two minima in the 
periodic state. At early times, the system's behavior is 
similar to the symmetric case, but at some time, spon- 
taneous symmetry-breaking occurs. The length of time 
before steady-state is reached in the asymmetric case de- 
pends on degree of noise introduced in the strength of 
fluid-substrate interaction 0. The fact that one of the 
peaks becomes weaker than the other indicates that in- 
stead of the arms growing equally from both substrates 
and simultaneously forming a bubble in the middle, the 
top arm, for example, grows faster and contributes more 
to the bubble (higher maximum) than the bottom arm 
(smaller maximum). But for the next bubble, the situ- 
ation is reversed, so the whole period encompasses both 
maxima; this period is roughly twice that of the symmet- 
ric case. 

We also examined the system response as we change 
H dynamically. For each value of H , the steady-state 
value of Bi, or the maximum and minimum in the os- 
cillatory regimes, are shown on the bifurcation diagram 
in Fig. 5. By abruptly increasing H from below Hi to 
above H§, the system switches from the symmetric time- 
independent low-velocity regime (as in Fig. 3a) to the 
time-independent high- velocity regime (as in Fig. 3e). 
Correspondingly, by abruptly decreasing H from above 
Hq to below Hi, the system switches from the Fig. 3e 
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to the Fig. 3a regime. But the transition region between 
these two points is highly complicated, and incorporates 
complex bifurcations between all the different states. For 
H\ < H < Hq, the behavior of the system depends on 
the starting value of H , the direction of change (increas- 
ing or decreasing) and the magnitude of the change. For 
example, as we gradually increase H from below Hi (see 
red curve in Fig 5), the time- independent symmetric be- 
havior (as in Fig. 3a) becomes unstable and symmetry 
breaking occurs at the bifurcation point H\, giving rise 
to the one of the two possible time-independent asym- 
metric types of patterns (as in Fig. 3b) . Once the sym- 
metry of the system is broken, we find that it is no longer 
possible to form arms that grow symmetrically from the 
top and bottom substrate and subsequently form bubbles 
that are symmetric with respect to top/bottom. Further 
gradual increases in H lead to the asymmetric oscilla- 
tions (as in Fig. 3d) for H > H 5 . For H 5 < H < H 6 , we 
only observe asymmetric oscillations, independent of the 
direction of change in H. Inside the parameter region 
H2 < H < H$, we can observe all the different types of 
behavior shown in Figs. 3b-d. Gradually decreasing H 
leads to a transition from asymmetric to symmetric os- 
cillations, and then to the time independent asymmetric 
state (blue curve) Figure 5 clearly shows that the 

system displays hystcrisis. 



Oscillations (both types) 
Time-independent 
Symmetric oscillations 
Asymmetric oscillations 
Time-independent 
Asymmetric oscillations 
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FIG. 5: Bifurcation diagram. Open circles mark asymmet- 
ric oscillations, filled circles represent symmetric oscillations, 
and straight lines indicate time-independent behavior (except 
green lines that depict both symmetric (H < H4) and asym- 
metric (H > H4) oscillations). All other parameters (except 
H) are the same as in Fig 2. In the oscillatory regimes, the 
top and bottom curves correspond to the maximum and min- 
imum of Bi. 



The non-decaying, time-periodic behavior in a simple 
binary fluid driven through the microchannel arises from 
interactions between the fluid and the patterned sub- 
strate, which introduces non-linearity into the dynami- 
cal system. These interactions act as an "activator" in 



reaction-diffusion systems |3J and are responsible for the 
positive feedback. We note that the periods of the os- 
cillations and the positions of the bifurcation points are 
dependent on the strength of the fluid-substrate interac- 
tions and the patch length. Thus, the system dynamics 
can potentially be controlled by varying these chemical 
features. In particular, other choices of patterns can po- 
tentially lead to new spatiotemporal patterns and dy- 
namical behavior. 
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